# Joint distribution F(x_1, x_2, ..., x_d) where x_k is the execution time of the k'th function


EmpiricalCDF <- NULL  
for (i in 1:nb.fct) {  
	EmpiricalCDF <- c(EmpiricalCDF, ecdf(mydata[,i]))  
}

ECDF <- function(i, x) EmpiricalCDF[[i]](x)

InvECDF <- function(j, p) {
	minimum <- -1
	for (i in 1:nb.run) {
		if (ECDF(j, mydata[i,j]) >= p) {
			if (mydata[i,j] < minimum || minimum == -1) minimum <- mydata[i,j]
		}
	}
	minimum
}

JointCDF <- function(v) {
	print(v)
	values <- NULL
	for (j in 1:nb.fct) {
		values <- c(values, ECDF(j, v[j]))
	}
	values <- as.matrix(values)
	values <- t(values)
	print(values)
	C.n(u=values, U=pobs(mydata)) # return N/nb.run where N is the number of runs for which the execution time of EVERY function X_i is <= v[i]
}

#mydata <- read.table("output_10_100.txt", header=TRUE, sep=" ")
mydata <- read.table("full.txt", header=TRUE, sep=" ")

nb.run <- dim(mydata)[1]
nb.fct <- dim(mydata)[2]


#test case
nb.run <- 30000
nb.fct <- 2
mydata = matrix(runif(nb.run * nb.fct), c(nb.run,nb.fct))
mydata <- mydata * 50
#end of test case

print(paste("The data:", nb.fct, "functions over",  nb.run, "runs", sep=" "))
#print(mydata)

N <- 2^nb.fct
I <- array(0, c(N, nb.fct))
sumI <- array(0, c(N))

for (i in 2:N) {
	I[i,] <- I[i-1,]
	j <- 1
	I[i, j] <- I[i, j] + 1
	while (I[i, j] > 1) {
		I[i, j] <- 0
		j <- j + 1
		I[i, j] <- I[i, j] + 1
	}
	sumI[i] <- sum(I[i,])
}

print(paste("N =", N, "= 2 ^ the number of functions/variables", sep=" "))
print(paste("The binary matrix I that counts from 1 to N, i.e. from 1 to", N, sep=" "))
print(I)

Hmultiplier <- array(0, c(N))
for (i in 1:N) {
	if ((nb.fct  -  sumI[i]) %% 2 == 0) Hmultiplier[i] <- 1
	else  Hmultiplier[i] <- -1
}
	
print("The vector of multipliers used in the computation of the HyperCube, in {-1, 1}")
print(Hmultiplier)


HyperCube <- function (b, h) {
	if (length(b) != nb.fct) {
		print("Error in Function HyperCube: argument too long")
	}
	samplePoints <- array(0, c(N, nb.fct))
	Uval <- array(0, c(N, nb.fct))
	for (i in 1:N) {
		samplePoints[i,] <- b + (h * I[i,])
		for (j in 1:nb.fct) Uval[i,j] <- ECDF(j, samplePoints[i,j])	
	}
	#print(samplePoints)
	#print(Uval)
	H <- C.n(u=Uval, U=pobs(mydata))	
	#print("H = ")
	#print(H)
	#print("Hmultiplier = ")
	#print(Hmultiplier)
	Hm <- Hmultiplier * H
	#print("Hm = ")
	#print(Hm)
	#print(sum(Hm))
	sum(Hm)
}

x <- c(0,0)
y <- HyperCube(x, 25)

#observed.WCET <- rowSums(mydata)
#observed.WCET.ecdf <- ecdf(observed.WCET)
#plot(observed.WCET.ecdf)

#sample.points.U <- as.matrix(expand.grid(rep(list(0:5), nb.fct)))
#sample.points.U <- (sample.points.U / 5)
#nb.points <- nrow(sample.points.U)

#sample.points.X <- array(0, c(nb.points, nb.fct))
#for (i in 1:nb.points) {
#	for (j in 1:nb.fct) {
#		sample.points.X[i,j] <- InvECDF(j, sample.points.U[i,j])
#	}
#}
#sample.WCET <- rowSums(sample.points.X)
#estimated.ecdf <- C.n(u=sample.points.U, U=pobs(mydata))


#sample.points <- matrix(runif(nb.point * nb.fct), nb.point, nb.fct)
#estimated.copula.ecdf <- C.n(u=sample.points, U=pobs(mydata))
 
#y <- dCn(u=mat, U=pobs(mydata), j.ind=1:nb.fct, b=1/sqrt(nrow(pobs(mydata))))